Create an interaction variable for time_point and sample
merged_data\(interaction_var <- interaction(merged_data\)sample, merged_data$time_point)
ggplot(merged_data, aes(x=interaction_var,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
geom_bar(stat=“identity”, colour=“white”, linewidth=0.1) + #plot stacked bars with white borders
scale_fill_manual(values=phylum_colors) +
facet_nested(. ~ type, scales=“free”) + #facet per day and treatment
guides(fill = guide_legend(ncol = 1)) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
axis.title.x = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(linewidth = 0.5, linetype = “solid”, colour = “black”)) +
labs(fill=“Phylum”,y = “Relative abundance”,x=“Samples”)+
scale_x_discrete(labels = function(x) gsub(“.\.”, ””, x)) +
facet_wrap(~type, scales = ”free”, labeller = as_labeller(function(label) gsub(”.\.”, ““, label))) #only show time_point label
### Phylum relative abundances
```{.r .script-source}
phylum_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
left_join(genome_metadata, by = join_by(genome == genome)) %>%
group_by(sample,phylum) %>%
summarise(relabun=sum(count))
phylum_summary %>%
group_by(phylum) %>%
summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
arrange(-mean) %>%
tt()
tinytable_hskwrb5m3uj3l6n4y8i5
| phylum |
mean |
sd |
| p__Bacteroidota |
0.376352935 |
0.205078492 |
| p__Bacillota_A |
0.250845368 |
0.165219604 |
| p__Bacillota |
0.116698680 |
0.146369286 |
| p__Pseudomonadota |
0.095767942 |
0.159976047 |
| p__Campylobacterota |
0.052429876 |
0.092626270 |
| p__Verrucomicrobiota |
0.029518627 |
0.072768229 |
| p__Desulfobacterota |
0.023304680 |
0.036356268 |
| p__Actinomycetota |
0.012917195 |
0.107880476 |
| p__Chlamydiota |
0.010449046 |
0.059349062 |
| p__Fusobacteriota |
0.010359534 |
0.028167393 |
| p__Cyanobacteriota |
0.009098787 |
0.016424861 |
| p__Bacillota_C |
0.004658039 |
0.006626014 |
| p__Spirochaetota |
0.003962784 |
0.012243051 |
| p__Bacillota_B |
0.002436629 |
0.004905803 |
| p__Elusimicrobiota |
0.001199879 |
0.006464719 |
phylum_arrange <- phylum_summary %>%
group_by(phylum) %>%
summarise(mean=mean(relabun)) %>%
arrange(-mean) %>%
select(phylum) %>%
pull()
phylum_summary %>%
filter(phylum %in% phylum_arrange) %>%
mutate(phylum=factor(phylum,levels=rev(phylum_arrange))) %>%
ggplot(aes(x=relabun, y=phylum, group=phylum, color=phylum)) +
scale_color_manual(values=phylum_colors[rev(phylum_arrange)]) +
geom_jitter(alpha=0.5) +
theme_minimal() +
theme(legend.position="none") +
labs(y="Phylum",x="Relative abundance")

Taxonomy boxplot
Family
family_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
group_by(sample,family) %>%
summarise(relabun=sum(count))
family_summary %>%
group_by(family) %>%
summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
arrange(-mean) %>%
tt()
tinytable_xlfq8zdtmm6vql0bajvc
| family |
mean |
sd |
| f__Bacteroidaceae |
2.189238e-01 |
0.1404298747 |
| f__Lachnospiraceae |
1.389745e-01 |
0.1089608538 |
| f__Tannerellaceae |
1.016053e-01 |
0.0802366634 |
| f__Helicobacteraceae |
5.197482e-02 |
0.0922005028 |
| f__Mycoplasmoidaceae |
3.651280e-02 |
0.0753015455 |
| f__Erysipelotrichaceae |
3.461428e-02 |
0.0450141949 |
| f__UBA3700 |
3.378612e-02 |
0.0555421417 |
| f__Marinifilaceae |
2.741603e-02 |
0.0272043836 |
| f__ |
2.740335e-02 |
0.0833966811 |
| f__DTU072 |
2.668734e-02 |
0.0526713763 |
| f__Rikenellaceae |
2.665357e-02 |
0.0462652505 |
| f__Enterobacteriaceae |
2.660196e-02 |
0.0913429878 |
| f__Coprobacillaceae |
2.521928e-02 |
0.0887342893 |
| f__Desulfovibrionaceae |
2.330468e-02 |
0.0363562684 |
| f__Ruminococcaceae |
1.838184e-02 |
0.0427420605 |
| f__LL51 |
1.738859e-02 |
0.0683714736 |
| f__Rhizobiaceae |
1.578130e-02 |
0.0766280394 |
| f__UBA3830 |
1.560351e-02 |
0.0451659123 |
| f__Mycobacteriaceae |
1.228390e-02 |
0.1079320546 |
| f__Akkermansiaceae |
1.213003e-02 |
0.0315247341 |
| f__Chlamydiaceae |
1.044905e-02 |
0.0593490616 |
| f__Fusobacteriaceae |
1.035953e-02 |
0.0281673926 |
| f__CAG-239 |
9.034415e-03 |
0.0150424841 |
| f__Enterococcaceae |
8.322037e-03 |
0.0463917038 |
| f__Gastranaerophilaceae |
7.638375e-03 |
0.0143794214 |
| f__Oscillospiraceae |
6.565519e-03 |
0.0077974910 |
| f__UBA1997 |
6.303806e-03 |
0.0307918484 |
| f__Streptococcaceae |
6.291980e-03 |
0.0340225424 |
| f__UBA1242 |
4.618700e-03 |
0.0152906152 |
| f__Brevinemataceae |
3.962784e-03 |
0.0122430505 |
| f__Acutalibacteraceae |
3.335081e-03 |
0.0108974502 |
| f__UBA660 |
3.159125e-03 |
0.0138778160 |
| f__Clostridiaceae |
2.808963e-03 |
0.0170355922 |
| f__RUG11792 |
2.784774e-03 |
0.0249664653 |
| f__Peptococcaceae |
2.436629e-03 |
0.0049058034 |
| f__MGBC116941 |
2.134904e-03 |
0.0093302936 |
| f__Acidaminococcaceae |
1.921070e-03 |
0.0050068934 |
| f__CAG-508 |
1.786832e-03 |
0.0063826540 |
| f__Anaerovoracaceae |
1.551174e-03 |
0.0036295746 |
| f__Moraxellaceae |
1.468042e-03 |
0.0096884882 |
| f__RUG14156 |
1.460412e-03 |
0.0045605184 |
| f__Staphylococcaceae |
1.345783e-03 |
0.0050584305 |
| f__Elusimicrobiaceae |
1.199879e-03 |
0.0064647185 |
| f__CAG-288 |
9.379861e-04 |
0.0059852367 |
| f__Anaerotignaceae |
8.884602e-04 |
0.0040242420 |
| f__CALVMC01 |
7.428930e-04 |
0.0043360033 |
| f__Eggerthellaceae |
6.332937e-04 |
0.0021152288 |
| f__Massilibacillaceae |
6.162148e-04 |
0.0016267915 |
| f__UBA1820 |
4.680638e-04 |
0.0012990360 |
| f__Arcobacteraceae |
4.550604e-04 |
0.0050029784 |
| f__CAG-274 |
4.466884e-04 |
0.0021903937 |
| f__Burkholderiaceae_C |
3.656162e-04 |
0.0047810524 |
| f__Muribaculaceae |
3.575580e-04 |
0.0009726557 |
| f__UBA932 |
3.379555e-04 |
0.0011520742 |
| f__Hepatoplasmataceae |
2.954146e-04 |
0.0038630477 |
| f__Rhodobacteraceae |
2.924483e-04 |
0.0038242573 |
| f__Weeksellaceae |
2.739210e-04 |
0.0031292133 |
| f__Eubacteriaceae |
1.627561e-04 |
0.0006691724 |
| f__Sphingobacteriaceae |
1.488164e-04 |
0.0012387571 |
| f__Devosiaceae |
1.472568e-04 |
0.0015006126 |
| f__Pumilibacteraceae |
1.262477e-04 |
0.0007602888 |
| f__WRAU01 |
9.491039e-05 |
0.0012411144 |
| f__Peptostreptococcaceae |
2.260586e-05 |
0.0002956099 |
family_arrange <- family_summary %>%
group_by(family) %>%
summarise(mean=sum(relabun)) %>%
arrange(-mean) %>%
select(family) %>%
pull()
family_summary %>%
left_join(genome_metadata %>% select(family,phylum) %>% unique(),by=join_by(family==family)) %>%
left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
filter(family %in% family_arrange[1:20]) %>%
mutate(family=factor(family,levels=rev(family_arrange[1:20]))) %>%
filter(relabun > 0) %>%
ggplot(aes(x=relabun, y=family, group=family, color=phylum)) +
scale_color_manual(values=phylum_colors[-8]) +
geom_jitter(alpha=0.5) +
facet_grid(.~type)+
theme_minimal() +
labs(y="Family", x="Relative abundance", color="Phylum")

Genus
genus_summary <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
left_join(genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
group_by(sample,genus) %>%
summarise(relabun=sum(count)) %>%
filter(genus != "g__")
genus_summary %>%
group_by(genus) %>%
summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
arrange(-mean) %>%
tt()
tinytable_rr62kpha5pxiisbuthf1
| genus |
mean |
sd |
| g__Bacteroides |
1.337013e-01 |
0.0929327655 |
| g__Parabacteroides |
9.568074e-02 |
0.0803866834 |
| g__Phocaeicola |
6.854013e-02 |
0.0793558268 |
| g__Mycoplasmoides |
3.015602e-02 |
0.0749534641 |
| g__Helicobacter_J |
2.973811e-02 |
0.0592188466 |
| g__Odoribacter |
2.522456e-02 |
0.0267412147 |
| g__Roseburia |
2.356827e-02 |
0.0556733585 |
| g__NHYM01 |
2.223671e-02 |
0.0792620634 |
| g__Alistipes |
2.182523e-02 |
0.0282922451 |
| g__Coprobacillus |
1.990305e-02 |
0.0873936500 |
| g__Agrobacterium |
1.578130e-02 |
0.0766280394 |
| g__Corynebacterium |
1.228390e-02 |
0.1079320546 |
| g__Akkermansia |
1.213003e-02 |
0.0315247341 |
| g__Fusobacterium_A |
1.026752e-02 |
0.0281722895 |
| g__Kineothrix |
8.698172e-03 |
0.0406735379 |
| g__Proteus |
8.556721e-03 |
0.0677872967 |
| g__Dielma |
8.477779e-03 |
0.0090773537 |
| g__CAG-95 |
8.014504e-03 |
0.0203997532 |
| g__JAAYNV01 |
7.901098e-03 |
0.0194443338 |
| g__UBA866 |
7.175409e-03 |
0.0291873132 |
| g__Desulfovibrio |
6.936071e-03 |
0.0210368286 |
| g__Enterococcus |
6.919414e-03 |
0.0453969341 |
| g__Ureaplasma |
6.356785e-03 |
0.0137373135 |
| g__Lactococcus |
6.291980e-03 |
0.0340225424 |
| g__Lacrimispora |
5.993828e-03 |
0.0097829951 |
| g__Parabacteroides_B |
5.924566e-03 |
0.0099792665 |
| g__CALXRO01 |
5.698293e-03 |
0.0306767185 |
| g__Citrobacter |
5.650169e-03 |
0.0332631355 |
| g__NSJ-61 |
5.495080e-03 |
0.0197757124 |
| g__Breznakia |
5.440147e-03 |
0.0235692957 |
| g__Clostridium_AQ |
5.263895e-03 |
0.0121113284 |
| g__Salmonella |
5.073541e-03 |
0.0183551010 |
| g__Bilophila |
4.925550e-03 |
0.0088098210 |
| g__Hungatella_A |
4.755524e-03 |
0.0095127458 |
| g__MGBC136627 |
4.313746e-03 |
0.0162110012 |
| g__Escherichia |
4.139378e-03 |
0.0264569200 |
| g__UMGS1251 |
4.111189e-03 |
0.0072426766 |
| g__Hungatella |
4.068243e-03 |
0.0189714607 |
| g__Clostridium_Q |
3.965072e-03 |
0.0052000165 |
| g__Brevinema |
3.962784e-03 |
0.0122430505 |
| g__Thomasclavelia |
3.856935e-03 |
0.0108480346 |
| g__Scatousia |
3.607252e-03 |
0.0102197530 |
| g__Enterocloster |
3.578040e-03 |
0.0047221634 |
| g__Mailhella |
3.569833e-03 |
0.0101940649 |
| g__Copromonas |
3.557270e-03 |
0.0050035176 |
| g__Ventrimonas |
3.472264e-03 |
0.0070914148 |
| g__Caccovivens |
3.304213e-03 |
0.0121527504 |
| g__Lawsonia |
3.250547e-03 |
0.0116543597 |
| g__Fournierella |
3.179598e-03 |
0.0062039376 |
| g__Limenecus |
3.126281e-03 |
0.0065432982 |
| g__MGBC133411 |
3.006834e-03 |
0.0074209199 |
| g__Mucinivorans |
2.866176e-03 |
0.0371005376 |
| g__Sarcina |
2.808963e-03 |
0.0170355922 |
| g__Acetatifactor |
2.706113e-03 |
0.0058179611 |
| g__Eisenbergiella |
2.665864e-03 |
0.0067990748 |
| g__Bacteroides_G |
2.651345e-03 |
0.0344120334 |
| g__CAJLXD01 |
2.603014e-03 |
0.0087025932 |
| g__Blautia |
2.528781e-03 |
0.0061107674 |
| g__C-19 |
2.248901e-03 |
0.0048466847 |
| g__Butyricimonas |
2.191467e-03 |
0.0049842986 |
| g__Velocimicrobium |
2.175479e-03 |
0.0066412559 |
| g__CAZU01 |
2.170696e-03 |
0.0065598332 |
| g__MGBC116941 |
2.134904e-03 |
0.0093302936 |
| g__Intestinimonas |
2.057548e-03 |
0.0039214573 |
| g__Negativibacillus |
2.044877e-03 |
0.0054857591 |
| g__Rikenella |
1.962168e-03 |
0.0036910712 |
| g__Phascolarctobacterium |
1.921070e-03 |
0.0050068934 |
| g__RGIG6463 |
1.768909e-03 |
0.0039495868 |
| g__JALFVM01 |
1.721982e-03 |
0.0038441918 |
| g__Oscillibacter |
1.492792e-03 |
0.0024898538 |
| g__Pseudoflavonifractor |
1.485397e-03 |
0.0027590442 |
| g__Acinetobacter |
1.468042e-03 |
0.0096884882 |
| g__Citrobacter_A |
1.376632e-03 |
0.0060011629 |
| g__Staphylococcus |
1.345783e-03 |
0.0050584305 |
| g__RGIG4733 |
1.288541e-03 |
0.0040266511 |
| g__UBA1436 |
1.199879e-03 |
0.0064647185 |
| g__Lachnotalea |
1.194066e-03 |
0.0048896855 |
| g__Ruthenibacterium |
1.187962e-03 |
0.0053331515 |
| g__14-2 |
1.170787e-03 |
0.0095739735 |
| g__Beduini |
1.160220e-03 |
0.0025025910 |
| g__Scatocola |
1.108079e-03 |
0.0044726561 |
| g__Faecisoma |
1.072888e-03 |
0.0055281036 |
| g__Enterococcus_A |
1.071313e-03 |
0.0098572488 |
| g__Faecimonas |
9.764964e-04 |
0.0082595968 |
| g__RGIG9287 |
9.526092e-04 |
0.0092684203 |
| g__CAG-345 |
9.379861e-04 |
0.0059852367 |
| g__Blautia_A |
9.100332e-04 |
0.0028927770 |
| g__RGIG1896 |
8.246024e-04 |
0.0062408705 |
| g__CAG-269 |
7.889086e-04 |
0.0046906408 |
| g__Marseille-P3106 |
7.845423e-04 |
0.0017535790 |
| g__WRHT01 |
6.354364e-04 |
0.0026829604 |
| g__Eggerthella |
6.332937e-04 |
0.0021152288 |
| g__IOR16 |
6.324101e-04 |
0.0020540977 |
| g__Anaerotruncus |
6.192165e-04 |
0.0016862064 |
| g__RUG14156 |
6.079829e-04 |
0.0022026464 |
| g__CHH4-2 |
6.073170e-04 |
0.0019890673 |
| g__Serratia_A |
5.792071e-04 |
0.0075741155 |
| g__CAG-56 |
4.857941e-04 |
0.0016254203 |
| g__Merdimorpha |
4.680638e-04 |
0.0012990360 |
| g__MGBC140009 |
4.624605e-04 |
0.0023930648 |
| g__CALURL01 |
4.581073e-04 |
0.0016646244 |
| g__Aliarcobacter |
4.550604e-04 |
0.0050029784 |
| g__Scatenecus |
4.467347e-04 |
0.0019650337 |
| g__RGIG8482 |
4.347613e-04 |
0.0029582243 |
| g__Enterobacter |
4.025794e-04 |
0.0041076314 |
| g__Klebsiella |
4.007019e-04 |
0.0048624260 |
| g__Caccenecus |
3.895102e-04 |
0.0017702443 |
| g__Alcaligenes |
3.656162e-04 |
0.0047810524 |
| g__Plesiomonas |
3.590755e-04 |
0.0026947991 |
| g__HGM05232 |
3.575580e-04 |
0.0009726557 |
| g__JAHHSE01 |
3.550107e-04 |
0.0014818046 |
| g__Egerieousia |
3.379555e-04 |
0.0011520742 |
| g__Emergencia |
3.373704e-04 |
0.0017351293 |
| g__Enterococcus_B |
3.313107e-04 |
0.0022138492 |
| g__Stoquefichus |
2.990680e-04 |
0.0020385612 |
| g__Hepatoplasma |
2.954146e-04 |
0.0038630477 |
| g__Paracoccus |
2.924483e-04 |
0.0038242573 |
| g__Moheibacter |
2.739210e-04 |
0.0031292133 |
| g__Scatomorpha |
2.610126e-04 |
0.0010128258 |
| g__UBA7185 |
2.405856e-04 |
0.0014474682 |
| g__Eubacterium |
1.627561e-04 |
0.0006691724 |
| g__Sphingobacterium |
1.488164e-04 |
0.0012387571 |
| g__Devosia |
1.472568e-04 |
0.0015006126 |
| g__Anaerosporobacter |
1.437105e-04 |
0.0012673026 |
| g__Caccomorpha |
1.366945e-04 |
0.0010479478 |
| g__UBA2658 |
1.292159e-04 |
0.0007163845 |
| g__Protoclostridium |
1.262477e-04 |
0.0007602888 |
| g__Angelakisella |
1.253643e-04 |
0.0009168117 |
| g__Cetobacterium_A |
9.201325e-05 |
0.0008667719 |
| g__Rahnella |
6.395024e-05 |
0.0008362579 |
| g__Peptostreptococcus |
2.260586e-05 |
0.0002956099 |
genus_arrange <- genus_summary %>%
group_by(genus) %>%
summarise(mean=sum(relabun)) %>%
filter(genus != "g__")%>%
arrange(-mean) %>%
select(genus) %>%
mutate(genus= sub("^g__", "", genus)) %>%
pull()
genus_summary %>%
left_join(genome_metadata %>% select(genus,phylum) %>% unique(),by=join_by(genus==genus)) %>%
left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
mutate(genus= sub("^g__", "", genus)) %>%
filter(genus %in% genus_arrange[1:20]) %>%
mutate(genus=factor(genus,levels=rev(genus_arrange[1:20]))) %>%
filter(relabun > 0) %>%
ggplot(aes(x=relabun, y=genus, group=genus, color=phylum)) +
scale_color_manual(values=phylum_colors[-c(3,4,6,8)]) +
geom_jitter(alpha=0.5) +
facet_grid(.~type)+
theme_minimal() +
labs(y="Family", x="Relative abundance", color="Phylum")

Alpha diversity
# Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 0) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(richness = 1) %>%
rownames_to_column(var = "sample")
neutral <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(neutral = 1) %>%
rownames_to_column(var = "sample")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, tree = genome_tree) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(phylogenetic = 1) %>%
rownames_to_column(var = "sample")
# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
to.elements(., GIFT_db) %>%
traits2dist(., method = "gower")
functional <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, dist = dist) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(functional = 1) %>%
rownames_to_column(var = "sample") %>%
mutate(functional = if_else(is.nan(functional), 1, functional))
# Merge all metrics
alpha_div <- richness %>%
full_join(neutral, by = join_by(sample == sample)) %>%
full_join(phylogenetic, by = join_by(sample == sample)) %>%
full_join(functional, by = join_by(sample == sample))
Wild samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="0_Wild") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b", "#6b7398")) +
scale_fill_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b50", "#6b739850")) +
facet_wrap(. ~ metric ) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)

Acclimation samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="1_Acclimation") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b", "#6b7398")) +
scale_fill_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b50", "#6b739850")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)

Antibiotics samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="2_Antibiotics") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b", "#6b7398")) +
scale_fill_manual(name="species",
breaks=c("Podarcis_liolepis","Podarcis_muralis"),
labels=c("Podarcis liolepis","Podarcis muralis"),
values=c("#e5bd5b50", "#6b739850")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)

Transplant_1 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="3_Transplant1") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)

Transplant_2 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="4_Transplant2") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)

Post-Transplant_1 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="5_Post-FMT1") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)

Post-Transplant_2 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="6_Post-FMT2") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b183","#d57d2c", "#4477AA")) +
scale_fill_manual(name="type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Control","Hot control", "Treatment"),
values=c("#76b18350","#d57d2c50", "#4477AA50")) +
facet_wrap(. ~ metric) +
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
)

Beta diversity
beta_q0n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 0)
beta_q1n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1)
beta_q1p <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1, tree = genome_tree)
beta_q1f <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1, dist = dist)
Permanovas
Load required data
meta <- column_to_rownames(sample_metadata, "Tube_code")
1. Are the wild populations similar?
Wild: P.muralis vs P.liolepis
[1] TRUE
Number of samples used
[1] 28
beta_div_richness_wild<-hillpair(data=wild.counts, q=0)
beta_div_neutral_wild<-hillpair(data=wild.counts, q=1)
beta_div_phylo_wild<-hillpair(data=wild.counts, q=1, tree=genome_tree)
beta_div_func_wild<-hillpair(data=wild.counts_func, q=1, dist=dist)
Richness
| species |
1 |
1.631953 |
0.207198 |
6.795072 |
0.001 |
| Residual |
26 |
6.244347 |
0.792802 |
NA |
NA |
| Total |
27 |
7.876300 |
1.000000 |
NA |
NA |
Neutral
| species |
1 |
2.018797 |
0.2566342 |
8.976049 |
0.001 |
| Residual |
26 |
5.847641 |
0.7433658 |
NA |
NA |
| Total |
27 |
7.866438 |
1.0000000 |
NA |
NA |
Phylogenetic
| species |
1 |
0.3786052 |
0.2108638 |
6.947419 |
0.001 |
| Residual |
26 |
1.4168908 |
0.7891362 |
NA |
NA |
| Total |
27 |
1.7954960 |
1.0000000 |
NA |
NA |
Functional
| species |
1 |
0.0787916 |
0.1594096 |
4.930642 |
0.062 |
| Residual |
26 |
0.4154800 |
0.8405904 |
NA |
NA |
| Total |
27 |
0.4942716 |
1.0000000 |
NA |
NA |

2. Do the antibiotics work?
Acclimation vs antibiotics
[1] TRUE
Number of samples used
[1] 55
beta_div_richness_treat<-hillpair(data=treat.counts, q=0)
beta_div_neutral_treat<-hillpair(data=treat.counts, q=1)
beta_div_phylo_treat<-hillpair(data=treat.counts, q=1, tree=genome_tree)
beta_div_func_treat<-hillpair(data=treat.counts_func, q=1, dist=dist)
Richness
| time_point |
1 |
2.042510 |
0.0911466 |
6.563244 |
0.001 |
| species |
1 |
2.300885 |
0.1026765 |
7.393486 |
0.001 |
| individual |
26 |
9.974365 |
0.4451038 |
1.232725 |
0.003 |
| Residual |
26 |
8.091314 |
0.3610731 |
NA |
NA |
| Total |
54 |
22.409075 |
1.0000000 |
NA |
NA |
Neutral
| time_point |
1 |
2.189257 |
0.1017464 |
8.979451 |
0.001 |
| species |
1 |
3.046568 |
0.1415902 |
12.495800 |
0.001 |
| individual |
26 |
9.941985 |
0.4620568 |
1.568386 |
0.001 |
| Residual |
26 |
6.338992 |
0.2946066 |
NA |
NA |
| Total |
54 |
21.516802 |
1.0000000 |
NA |
NA |
Phylogenetic
| time_point |
1 |
2.1642230 |
0.2000313 |
18.460841 |
0.001 |
| species |
1 |
0.8686306 |
0.0802844 |
7.409427 |
0.001 |
| individual |
26 |
4.7385071 |
0.4379630 |
1.554596 |
0.004 |
| Residual |
26 |
3.0480626 |
0.2817214 |
NA |
NA |
| Total |
54 |
10.8194233 |
1.0000000 |
NA |
NA |
Functional
| time_point |
1 |
2.2832796 |
0.3550091 |
33.1502018 |
0.001 |
| species |
1 |
0.0227598 |
0.0035387 |
0.3304425 |
0.712 |
| individual |
26 |
2.3347730 |
0.3630154 |
1.3037622 |
0.208 |
| Residual |
26 |
1.7907967 |
0.2784368 |
NA |
NA |
| Total |
54 |
6.4316092 |
1.0000000 |
NA |
NA |
beta_richness_nmds_treat <- beta_div_richness_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = c("sample" = "Tube_code"))
beta_neutral_nmds_treat <- beta_div_neutral_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = c("sample" = "Tube_code"))
beta_phylo_nmds_treat <- beta_div_phylo_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = join_by(sample == Tube_code))
beta_func_nmds_treat <- beta_div_func_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = join_by(sample == Tube_code))
